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ABSTRACT 

We introduce a powerful semi-numeric modeling tool, 21cmFAST, designed to efficiently sim- 
ulate the cosmological 21-cm signal. Our code generates 3D realizations of evolved density, 
ionization, peculiar velocity, and spin temperature fields, which it then combines to compute 
the 21-cm brightness temperature. Although the physical processes are treated with approxi- 
mate methods, we compare our results to a state-of-the-art large-scale hydrodynamic simula- 
tion, and find good agreement on scales pertinent to the upcoming observations (> 1 Mpc). 
The power spectra from 21cmFAST agree with those generated from the numerical simulation 
to within 10s of percent, down to the Nyquist frequency. We show results from a 1 Gpc sim- 
ulation which tracks the cosmic 21-cm signal down from z = 250, highlighting the various 
interesting epochs. Depending on the desired resolution, 21cmFAST can compute a redshift 
realization on a single processor in just a few minutes. Our code is fast, efficient, customizable 
and publicly available, making it a useful tool for 21-cm parameter studies. 

Key words: cosmology: theory - intergalactic medium - large scale structure of universe - 
early Universe - galaxies: formation - high-redshift - evolution 



1 INTRODUCTION 

Through challenging observational efforts, the high-redshift fron- 
tier has been incrementally pushed back in recent years. Glimpses 
of the z 6-8 Universe were provided by quasars (e.g. Ipan et al 
2006), candidate Lyman b reak galaxies ( e.g. iBouwens et al. 



2009) 



zUUd) . canaiaate Lyman b reak galaxies ( e.g. IBouwens 
20081 : iMcLure et all l2009l : iBouwens et '2009"; 'Ouch 



Lyman alpha em itters (LAEs; e.g. IShimasaku et alj 
Kashikawa et al. 2006), an d GRBs (e.g. [Cusumano et al 



etal. 




Greiner etaP 2009: S alvaterra et al.l l2009h . Unfortunately, these 



precious observations currently provide only a limited set of rel- 
atively bright objects. Luckily, we will soon be inundated with ob- 
servations probing this and even earlier epochs. These observations 
should include infrared spectra from the James Webb Space Tele- 
scope (JWST), the Thirty Meter Telescope (TMT), the Giant Mag- 
ellan Telescope (GMT), the European Extremely Large Telescope 
(E-ELT), wide-field LAE surveys from the Subaru HyperSupreme- 
Cam, as well as the E-mode CMB polarization power spectrum 
measured by the Planck satellite. Some of the most important infor- 
mation will come in the form of the redshifted 21-cm line of neu- 
tral hydrogen. Several interferometers will attempt to observe the 
cosmological 21-cm signal, including the Mileura Wide Field Ar- 
ray (MWA; .Bowman et al...2005.fl the Low Frequency Array (LO- 



FArA the Giant Metrewave Radio Telescope (GMRT; |Pen et all 
2008), the Precision A rray to Probe the Epoch of Reionization (PA- 
PER; P arsons et al E0091 and eventually the Square Kilometer Ar- 
ray (skaS The first generation of these instruments, most notably 
LOFAR and MWA, are not only scheduled to become operational 
within a year, but should also yield insight into the 3D distribution 
of HI, provided the systematics can be overcome (see th e recent 
reviews of Furlanetto et al. 200d : [Morales & Wvithdl2009h . 

However, interpreting this data will be quite challenging and 
no-doubt controversial, as foreshadowed by the confusion sur- 
rounding the scant, currently- available observations. There are two 
main challenges to overcome: (1) an extremely large parameter 
space, due to our poor understanding of the high-redshift Universe; 
(2) an enormous dynamic range (i.e. range of relevant scales). 

Theoretically, the dawn of the first astrophysical objects and 
reionization could be modeled from first principles using numerical 
simulations, which include the complex interplay of many physical 
processes. In practice however, simulating these epochs requires 
enormous simulation boxes. Gigaparsec scales are necessary to sta- 
tistically model ionized regions and absorption systems or create 
accurate mock spectra from the very rare high-redshift quasars. 
However, the simulations also require high enough resolution to re- 
solve the underlying sources and sinks of ionizing photons and the 
complex small-scale feedback mechanisms which regulate them. 
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Thus one is forced to make compromises: deciding which physical 
processes can be ignored, and how the others can be parameterized 
and efficiently folded into large-scale models. Furthermore, large- 
scale simulations are computationally costly (even when they sacri- 
fice completeness for speed by ignoring hydrodynamic processes) 
and thus are inefficient in large parameter studies. 

On the other hand, analytic models, while more approximate, 
are fast and can provide physical insight into the import of vari- 
ous processes. However, analytical models are hard-pressed to go 
beyond the linear regime, and beyond making fairly simple predic- 
tions such as the mean 21 -cm signal Ipurlanetto 2006), the prob- 
ability densit y function (PDF ; Furlanetto et al. 2004a) and power 
spectrum ( Pr itchard & Furlanetto 2007 ; Barkana 2009). The 21- 
cm tomographic signal should be rich in information, accommo- 
dating many additional, higher-order statistical probes, such as the 
bi- spectrum (Pritchard et al., in preparation), the difference PDF 
tearkana & Loebll2008h . etc. 

In this paper, we follow a path of compromise, attempting 
to preserve the most useful elements of both analytic and nu- 
meric approaches. We introduce a self-consistent, semi-numericafl 
simulation, specifically optimized to predict the high-redshift 21- 
cm signal. Through a combination of the excursion- set formal- 
ism and perturbation theory, our code can generate full 3D re- 
alizations of the density, ionization, velocity, spin temperature, 
and ultimately 21 -cm brightness temperature fields. Although the 
physical processes are treated with approximate methods, our 
results agree well with a state-of-the-art hydrodynamic simula- 
tion of reionization. However, unlike numerical simulations, re- 
alizations are computationally cheap and can be generated in a 
matter of minutes on a single processor, with modest memory 
requirements. Most importantly, our code is publicly available 
at http://www.astro.princeton .edu/^mesinger/Sim.html| We name 
our simulation 21cmFAST. 

Semi-Numerical a pproaches have al r eady proved invaluable 



in re i onization studies dZat 



20071: iGeil & Wvithd l2008l : 



net al. 2005; 



Mesinger & Fu rlanetto 
2009 : Choudhurv et al 



J Alvarez et al 

2OO9I : iThomas et al.ll2009h . Indeed, 21cmF AST is a more special 



ized version of our previous code, DexM ^ Mesinger & Furlanettol 
I2OO7I : hereafter MF07). The difference between the two is that 
21cmFAST bypasses the halo finding algorithm, resulting in a 
faster code with a larger dynamic range and more modest memory 
requirements. In this work, we also introduce some new additions 
to our code, mainly to compute the spin temperature. 

In ^ we compare predictions from 21cmFAST with those 
from hydrodynamic simulations of the various physical compo- 
nents comprising the 21 -cm signal in the post heating regime. Den- 
sity, ionization, peculiar velocity gradient, and full 21 -cm bright- 
ness temperature fields are explored in ^TTl ^2l2l ^231 ^2A\ re- 
spectively. In ^ we introduce our method for computing the spin 
temperature fields, with results from the complete calculation (in- 
cluding the spin temperature) presented in ^3.31 Finally in ^ we 
summarize our findings. 

Unless stated otherwise, we quote all quantities in comov- 
ing units. We adopt the background cosmological parameters (Qa, 
Qm, ^b, n, as. Ho) = (0.72, 0.28, 0.046, 0.96, 0.82, 70 km s"^ 
Mpc~^), matching t he five-year results of the WMAP satellite 
(iKomatsu et al.ll2009l) . 



^ By "semi-numerical" we mean using more approximate physics than nu- 
merical simulations, but capable of independently generating 3D realiza- 
tions. 



2 21-CM TEMPERATURE FLUCTUATIONS POST 
HEATING (Ts > T^) 

Our ultimate goal is to compute the 21 cm background, which re- 
quires a number of physics components. To identify them, note that 
the offset of the 21 -cm brightness temperature from the CMB tem- 
perature, Tj, along a line of sight (LOS) at observed frequency z/, 
can be written as fc.f. lFurlanetto et all2006l) : 



27XHl(l + ^nl 

1 + z 0.15 



dvr/dr + H J \ Ts 



10 Qmh^ 



1/2 



0.023 



mK, 



(1) 



where Ts is the gas spin temperature, t^q is the optical depth at the 
21-cm frequency z/q, ^ni(x, z) = p/p — lis the evolved (Eulerian) 
density contrast, H(z) is the Hubble parameter, dvr/dr is the co- 
moving gradient of the line of sight component of the comoving ve- 
locity, and all quantities are evaluated at redshift z = lyo/v — l. The 
final approximation makes the assumption that that dvr/dr <^ H, 
which is generally true for the pertinent redshifts and scales, though 
we shall return to this issue in ^2.31 

In this comparison section, we make the standard, simpli- 
fying assumption of working in the post-heating regime: Ts ^ 
Tj. For fiducial astrophysical models, this is likely a safe 
assumption during the bulk o f reionization jFurlanettol I2OO6I: 
Chen & Miralda-Escudd l2008l : ISantos et all l2008l : iBaek et all 



2009*). We will however revisit this assumption in ^ where we 
introduce our method for computing the spin temperature field. 

The remaining components of eq. [T]are the density, ^ni, the 
ionization, xhi, and the velocity gradient, dvr/dr. Below, we study 
these in turn, comparing 21 cmFAST to the hydrodynamic cosmo- 
logical simulation of Trac et al. (2008), using the same initial con- 
ditions (ICs). We perform "by-eye" comparisons at various red- 
shifts/stages of reionization, as well as one and two-point statis- 
tics: the PDFs (smoothed on several scales), and the power spectra. 
Since our code is designed to simulate the cosmological 21-cm sig- 
nal from neutral hydrogen, we study the regime before the likely 
completion of reionization, ^ > 7 (though pre sent data is even 
consistent with reionization completing at 2; <6; iLidz et al1l2007l : 
lMesingejl2009l) . ^ 

The simulations of iTrac et"aD (l2008h are the current "state-of- 
the-art" reionization simulations. They include simultaneous treat- 
ment of dark matter (DM) and gas, five-frequency radiative transfer 
(RT) on a 512^ grid, and they resolve Mhaio > 10^ M© ionizing 
sources with > 40 DM particles in a 143 Mpc box. 



2.1 Evolved Density Field 

We calculate the evolved density field in the same fashion as in the 
"parent" code, DexM, outlined in MF07. In short, we generate den- 
sity and velocity ICs in initial (Lagrangian) space, in roughly the 
same manner as numerical cosmological simulations. We then ap- 
proximate gravitational collapse by moving each initial matter par- 
ticle (whose mass is the total mass in the c orresponding IC c ell) ac- 
cording to first-order perturbation theory (IZerDovichlil970h . First- 
order perturbation theory is very computationally convenient as the 
displacement field is a separable function of space and time, so 
the spatial component need only be computed once for each re- 
alization/box. There is no separate treatment of baryons and DM. 
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Figure 1. Slices through the density field, A(x, z), for the gas, DM, linearly-extrapolated ICs, and 21cmFAST, at 2; = 7, (clockwise from top left). Each slice 
is 143 Mpc on a side and 0.19 Mpc thick. 



Readers interested in more details concerning this approach are en- 
couraged to check MF07. 

This a pproach to generating lar ge-s cale density fields w as also 
adopted bv IChoudhurv et all J2OO9I) and lSantos etaP ([2009), who 
briefly showed that the resulting fields at high- 2; traced the DM dis- 
tribution from an N-body code fairly well. Here we perform more 
extensive comparisons. The Lagrangian space ICs used for 21cm- 
FAST were initialized dX z — 300 on a 1536^ grid. The velocity 
fields used to perturb the ICs, as well as the resulting density fields 
presented below are 768^ We show results from both 768^ and 

^ Our code allows the ICs to be sampled onto a high-resolution, HI-RES - 
DIM^, grid, while the subsequent evolved density, ionization, etc. fields 
can be created at lower resolution, LOW-RES -DIM^. This allows efficient 
use of available memory, with the code storing at most HI-RES-DIM^ 
-h 4 X LOW-RES-DIM^ floating point numbers in RAM. However, the 



256^ boxes. We highlight here that it takes ^10 minutes to gener- 
ate the 768^ 21cmFAST density fleld from the 1536^ ICs at z = 7 
on a single processor Mac Pro desktop computer. To put this into 
perspective, a hydrodynamical simulation of this single realization 
would take approximately three days to run down to z = 7 on a 
1536-node supercomputing cluster. 

In Figure [TJ we show a slice through the evolved density fleld, 
A(x, = 1 + ^ni, at z = 7. Density flelds computed from the gas, 

Zel'Dovich perturbation approach, just as numerical simulations, requires 
that the evolved fields are adequately resolved by the high-resolution grid. 
Failure to do so can substantially underestimate the fluctuations in the den- 
sity fleld. We roughly flnd that the high-resolution grid should have cell 
sizes < 1 Mpc to accurately model density flelds at redshifts of interest 
{z < 40). For larger cell sizes, > 10 Mpc, the linear theory evolution 
option of 21cmFAST should be used. 
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Figure 2. PDFs of the density fields smoothed on scale i^fiiter^ computed from the gas (solid red curves), DM {dotted green curves), and 21cmFAST {dashed 
blue curves) fields, at z = 20, 15, 10, 7 {top to bottom). The left panel corresponds to i^fiiter = 0.5 Mpc; the right panel corresponds to i^fiiter = 5 Mpc. All 
smoothing was performed with a real-space, top-hat filter. 



DM, linearly-extrapolated ICs, and 21cmFAST are shown clock- 
wise from the top left panel. It is evident that the Zel'Dovich ap- 
proximation works quite well for this purpose, accurately repro- 
ducing the cosmic web. We do not capture baryonic physics, and 
so the 21cmFAST output looks more similar to the DM than the 
gas. However, hydrodynamics is not included even in most of the 



2006l : lMcOuinn 



present-day cosmological ( > 100 Mpc) simulations (e.g.llliev et al 
2006l : lMcOuinn et al.l200^ r 
200^ 



see the recent review in lXrac & Gnedin 



Also note that the linear density field is only accurate on large 
scales (> 10 Mpc at z ^ 7). Thus care should be taken in applying 
tools which rely on the linear density field (such as the standard ex- 
cursion set formalism) at smaller scales. Nevertheless, we include 
in 21cmFAST a feature to evolve the density field using fully lin- 
ear evolution, instead of the perturbation approach discussed above. 
This allows one to generate extremely large boxes at low resolution. 
When using this feature, one should make sure that the chosen cell 
size is indeed in the linear regime at the redshift of interest. Some 
results making use of this feature are presented below. 

In Fig. El we show the PDFs of the density fields smoothed on 
scale i?fiiter, computed from the gas (solid red curves), DM {dotted 
green curves), and 21cmFAST {dashed blue curves) fields, at 2; = 
20, 15, 10, 7 {top to bottom). From the left panel (i?fiiter = 0.5 
Mpc), we see that as structure formation progresses, we tend to 
increasingly over-predict the abundance of small scale underdensi- 
ties, and under-predict the abundance of large-scale overdensities. 
However, even dX z — 7, our PDFs are accurate at the percent level 
to over a dex around the mean density. Understandably, the agree- 
ment between the PDFs becomes better with increasing scale (see 
the right panel corresponding to i?fiiter = 5 Mpc). Interestingly, the 
DM distributions match the gas quite well, although this is some- 
what of a coincidence, as we shall see from the power spectra be- 
low. 

In Figure [S] we present the density power spectra, defined 
as A^5(A;,z) = /cV(27rV) (|^(k, z)!^)^. The solid red, dotted 
green, and dashed blue curves correspond to the gas, DM, and 
21cmFAST fields, respectively. On small scales (/c > 5 Mpc~^), 
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Figure 3. Density power spectra, A^^(k,z),of the gas {solid red curves), 
DM {dotted green curves), and 21cmFAST {dashed blue curves) fields, at 
at 2; = 20, 15, 10, 7 {top to bottom). 

the three fields have different power. The collapse of gas is ini- 
tially delayed with respect to the pressureless DM, resulting in less 
small scale power. The perturbation theory approach of 21cmFAST 
is closer in spirit to the DM evolution, but does not capture virial- 
ized structure. In fact the close agreement dX z — 7 between the 
gas and 21cmFAST density power spectra is a coincidence, with 
the small-scale flattening of the 21cmFAST power attributable to 
"shell-crossing" by the matter particles in the Zel'Dovich approxi- 
mation. During reionization, the evolution of the gas is very com- 
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plicated, since the power spectrum on small scales is sensitive to 
the th ermal history of the reionization model (e.g. iHui & GnedinI 

On large scales {k < 0.5 Mpc~^), all three power spectra 
agree remarkably at all epochs. To put this into perspective, neither 
the MWA nor LOFAR have sufficient sign al to noise to detec t the 
21-cm signal beyond /c > 2 Mpc~^ fe.g. lMcOuinn et al.ll2006h . 

2.2 Ionization Field 

We use a new, refi ned, semi-numeric algorithm, FFRT, presented in 
IZahn et al to generate ionization fields, xhi(x, z). The ion- 

ization fields h ave been exhaustiv ely compared against cosmologi- 
cal RT codes in lZahn etal yielding good agreement across 

a broad range of statistical diagnostics on moderate to large scales. 
Thus, we will not present further tests here. Instead, we merely out- 
line the procedure, and motivate some aspects with regards to the 
goals of 21cmFAST: speed and efficiency. 

We use t he excursion- set approach for identifying HII regions, 
pioneered bv ^ Furlanetto et al.l (12004b). The foundation of this ap- 
proach is to require that the number of ionizing photons inside a re- 
gion be larger than the number of hydrogen atoms it contains. Then 
ionized regions are identified via an excursion- set approach, start- 
ing at large scales and progressing to small scales, analogous to the 
derivation of the Press-Sch echter (PS) mass functions dBond et al] 
ll99ll : lLacev&Colelll993h . 

Specifically, we flag fully ionized cells in our box as those 
which meet the criteria /coii(x, i?) > C,~^ , where C is some 
efficiency parameter and /coii(x, i?) is the collapse fraction 
smoothed on decreasing scales, starting from a maximum i?max 
Mpc and going down to the cell size, Rce\\. Additionally, we al- 
low for partially-ionized cells by setting the cell's ionized frac- 
tion to C/coii(x, i?ceii) at the last filter step for those cells 
which are not fully ionizecQ. The ionizing photon horizon, i?max, 
is a free parameter which can be chosen to match the extrap- 
olated ionizing ph oton mean free path, in the ionized IGM, at 
z ^ 7-10 (e.g. Stor rie-Lombardi et al.ll99l : lMiralda-Escudel2003l : 
IChoudhurv et al.l 12008). The photon sinks dominating the mean 
free path of ionizing photons are likely too small to be resolved in 
reionization simulations. An effective horizon du e to photon sinks 
can delay the completion o f reionization (e.g. IChoudhurv et al.l 
l2009l : lFurlanetto & Mesinge r 2009), and cause a drop in large scale 
21-cm power, as we shall see below. 

There are two main differences between FFRT used in 21cm- 
FAST and the previous incarnation of our HII bubble finder used in 
DexM (MF07): (1) the use of the halo finder to generate ionization 
fields in DexM (MF07) vs. using the evolved density field and con- 
ditional PS to generate ionization fields in 21cmFAST; and (2) the 
bubble flagging algorithm, which in MF07 is taken to paint the en- 
tire spherical region enclosed by the filter as ionized ("flagging-the- 
entire- sphere"), whereas for 21cmFAST we just flag the central cell 
as ionized C'flag ging-the-central-celF'; for more information, see 
IZahn et al.lf200l Mesinger & Furlanettoll2007l and the appendix in 

^ Our algorithm also can optionally account for Poisson fluctuations 
in the halo number, when the mean collapse fraction becomes small, 
/coii(x, 2;, i^ceii) X ^ceii < SOMmin, whcrc Mceii is the total mass 
within the cell and our faintest ionizing sources correspond to a halo mass 
of Mniin- This last step is found to be somewhat important when the cell 
size increases to > 1 Mpc (see appendix in Zahn et al. 2010). This is left as 
an option since turning off such stochastic behavior allows the user to better 
track the deterministic redshift evolution of a single realization. 



IZahnetal.l ( [20TQh . These default settings of 21cmFAST were cho- 
sen to maximize speed and dynamic range, while minimizing the 
memory requirements. Nevertheless, they are left as user- adjustable 
options in the codes. 



The first difference noted above means that 21cmFAST does 
not explicitly resolve source halos. In MF07, we made use of 
a semi-numerically generated halo field, which accurately re- 
produces N-body halo fields down to non-linear scales (MF07; 
Mesinger et al., in preparation). As in numerical reionization sim- 
ulations, these halos were assumed to host ionizing sources. How- 
ever, the intermediate step of generating such halo source fields 
adds additional computation time, and generally requires many GB 
of RAM for typical cosmological uses. As for numerical simula- 
tions, this memory requirement means that simulation boxes are 
limited to < 200 Mpc if they wish to resolve atomically-cooled ha- 
los at z =7-10, and even smaller sizes if they wish to resolve these 
at higher redshifts or resolve molecularly-cooled halos. Although 
DexM's halo finder is much faster than N-body codes, and can gen- 
erate halo fields at a given redshift in a few hours on a single pro- 
cessor, extending the dynamic range even further without hundreds 
of GB of RAM would be very useful. Alternatives to extending 
the dynamic range have been proposed by McOuinn et al. ( 2007h : 
ISantos et all J2OO9I) . These involve stochastically populating cells 
with halos below the resolution threshold. Although computation- 
ally efficient, it is unclear if these alternatives preserve higher-order 
statistics of the non-Gaussian /coii(x, 2;, i?) field, as each cell is 
treated independently from the others. More fundamentally, the 
stochasticity involved makes it difficult to deterministically track 
the redshift evolution of a single realization: halos effectively pop 
in and out of existence from one redshift output to the next. 



Therefore, to increase speed and dynamic range, we use 
the FFRT algorithni , which uses the co nditional PS formalism 
(iLacev & Cold 1 19931 : ISomerville & Kolatti.1999) to generate the 
collapsed mass (i.e. ionizing source) field. Since conditional PS 
operates directly on the density field, without needing to resolve 
halos, one can have an enormous dynamic range with a relatively 
small loss in accuracy (compare FFRT and FFRT-S in IZahn et aP 
I2OIOI) . Most importantly, when computing /coii(x, 2;, i?) we use 
the non-linear density field, generated according to ^2.11 instead of 
the standard linear density field. The resulting ionization fields are 
a much better match to RT simulations than those generated from 
the linear density field dZahn et al.ll2oTQl : foreshadowed also by the 
ICs panel in Fig. [T] above). We normalize the resulting collapsed 
mass field to match the Sheth-Tormen (ST) mean collapse fraction, 
which in turn matches numerical simulations (see eq.[T4land asso- 
ciated discussion). 



The other major difference is that by default, FFRT in 21cm- 
FAST flags just the central filter cell, instead of the entire sphere 
enclosed by the filter, as in MF07. The main motivation for this 
switch is that the former algorithm is 0{N), while the later is 
slower: 0{N) at xhi ^ 1 but approaching 0{N'^) as xm 
. There are some other minor differences between the FFRT and 
the ionization scheme in MF07, such as the use of a sharp k- space 
filter instead of a spherical top-hat, but these have a smaller impact 
on the resulting ionization maps. 
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Figure 4. PDFs of the comoving LOS derivative of Vr [in units of H(z)], smoothed on scale i^fiiter = 0-5 Mpc (left) and 5.0 Mpc (right). Solid red curves 
are generated from the hydrodynamic simulation, while the dashed blue curves are generated by 21cmFAST. Redshifts corresponding to z =20, 15, 10, 7 
are shown top to bottom. All smoothing was performed with a real-space, top-hat filter. The dotted magenta curves were generated on comparable scales by 
21cmFAST with different initial conditions; however, they assume linear evolution of the density field, instead of the perturbation theory approach (see ^2. It . 



2.3 Peculiar Velocity Gradient Field 

Redshift space distortions, accounted for with the dvr/dr term in 
eq. are often ignored when simulating the 21 -cm signal. In 
the linear regime, redshift space distortions of the 21 -cm field are 
similar to the well-studied Kaiser effect (e.g. lKaise3ll987h . and the 
power spectrum of fluct uations is enhanced on a ll scales by a ge- 
ometri c factor of 1.87 (iBharadwai & Alii l2004l : iBarkana & Loebl 
l2005ah . However, the small scale overdensities, where redshift 
space distortions are most important, are also the regions whose 
21 -cm emission is first erased by "inside-out" reionization. Pre- 
liminary studies therefore concluded that redshift space distortions 
would only be noticea ble before reionization and in its early stages 
(lMcOuinnetalJl2006L MF07). As we are interested in accurately 
simulating the 21 -cm signal from all cosmological epochs, includ- 
ing pre-reionization, here we will compare the velocity gradient 
term from 21cmFAST and hydrodynamic simulations. 

Using the Zel'Dovich approximation on our 3D realizations, 
we can again efficiently move beyond the linear regime into the 
quasi-linear regime, and take into account correlations in the veloc- 
ity gradient field. In this first-order perturbation theory, the velocity 
field can be written as: 



v(k,z) 



2k 



(2) 



and so the derivative of the line-of- sight velocity, Vr where r for 
simplicity is oriented along a basis vector, can be written in k- space 
as: 



^ Note that this expression is exact, as long as the dvr/dr field is constant 
over the width of the 21 -cm line and dvr/dr <C H(z). Alternately, one 
can apply redshift space distortions when converting the comoving signal 
from the simulation box to an observed frequency. However, for sake of 
consistency, we perform all of our calculations in comoving space. 



dvr 
dr 



(k, Z) = ikrVri}^^ Z) 



I)(^)^nl(k), 



(3) 
(4) 



where differentiation is performed in k-space. The last approxima- 
tion is used for 21cmFAST, while the first, exact expression is used 
for the numerical simulatiorff]. 

In Fig.|4]we show the PDFs of the comoving LOS derivative 
of Vr [in units of if (z)], smoothed on scale i?fiiter =0.5 Mpc (left) 
and 5.0 Mpc (right). Solid red curves are generated from the hy- 
drodynamic simulation, while the dashed blue curves are generated 
by 21cmFAST. Redshifts corresponding to z =20, 15, 10, 7 are 
shown top to bottom. We see that our perturbation theory approach 
again does remarkably well in reproducing results from the hydro- 
dynamic simulation. The velocity gradients agree even better than 
the density fields, since the velocity field is coherent over larger 
scales. The shape of the distributions are noticeably non-linear on 
small scales and late times. The curves resemble PDFs of the sign- 
flipped non-linear density field, ^ni, which is understandable from 
eq. ©. The dotted magenta curves in the bottom right panels were 
generated on comparable scales by 21cmFAST with different ini- 
tial conditions; however, they assume linear evolution of the den- 
sity field, instead of the perturbation theory approach. As expected, 
linear evolution results in a symmetric Gaussian PDF. 

Do we reproduce the geometric, scale-free enhancement 
of the power spectrum on linear scales? In the top panel 
of Fig. O we plot dimensionless 21 -cm power spectra, 
Aii(A:,z) = A:V(27rV) (k, z)!^)^ where fei(x,z) = 
5Tb (x, z)/5Th{z) — 1. The spectra are generated by 21cmFAST in 



^ There is an inconsistency in the above equations for 21cmFAST, as eq. 
(|4) is applied to the non-linear density field, whereas eq. (|2) assumes a linear 
5. Nevertheless, as we shall show below, eq. (|4) reproduces the non-linear 
velocity gradient field from the numerical simulations remarkably well. 
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Figure 5. Top panel: Dimensionless power spectra from 21cmFAST, in a 
fully neutral universe, in the limit of Ts ^ Xy. The upper set of curves 
were computed including peculiar velocities, while the lower set were com- 
puted not including peculiar velocities. Bottom three panels: Ratios of the 
dimensionless power spectra with peculiar velocities to those not including 
peculiar velocities. The linear regime geometric enhancement of 1.87 is de- 
marcated by the upper horizontal dotted line. In all panels, the dot-dashed 
blue and dashed green curves were generated from the same ICs in a L = 1 
Gpc box, sampled with an initial resolution of Ax = 0.56 Mpc, with the 
evolved density, velocity, and ionization fields generated at lower resolu- 
tions of Ax = 3.3 Mpc (dashed green curves), and 10 Mpc (dot-dashed 
blue curves). The solid red curves were generated from a L = 5 Gpc box 
with a single resolution of Ax = 10 Mpc. However, the density field used 
for the solid red curves was generated assuming linear evolution, while the 
others were generated with first-order perturbation theory (see ^2.U . 

the limit of Ts ^ and assuming xhi = 1. The solid red curves 
correspond to a 5 Gpc box with Ax = 10 Mpc cells, while the dot- 
dashed blue and dashed green curves correspond to a 1 Gpc box 
with different resolutions. The upper set of curves were computed 
including peculiar velocities, while the lower set were computed 
not including peculiar velocities. The bottom three panels show the 
ratios of the power spectra that include redshift space distortions to 
those that do not. 

Indeed the red curves in Fig. [5j which were evolved linearly, 
accurately capture the enhancement factor of 1.87, shown with a 
dotted horizontal line. The other two curves, which include first 
order non-linear effects, show an enhancement of power in ex- 
cess of the purely geometric factor. From eq. ©, one sees that a 
high- value tail in the density distribution resulting from non-linear 
evolution would drive a corresponding negative tail in the dvr/dr 
distributions, which in turn enhances the 21 -cm signal through the 
{l/(dvr/dr/H) + 1) term in eq. [T] Although the dvr/dr distri- 
butions are zero-mean, the distributions of l/{dvr/dr/H + 1) are 
not. The bias to higher values is further enhanced when weighted 
by the local density as in the 6Tb expression, A/(dvr/dr/H + 1). 
Intuitively, infall in overdense regions causes photons emitted there 
to travel farther in order to reach a fixed relative redshift; there- 



Figure 6. Ratios of the dimensional power spectra, STt{z)^ A 21 {k, z), 
computed including peculiar velocities to those not including peculiar ve- 
locities. Panels correspond to (z,xm) = (9.00, 0.86), (7.73, 0.65), (7.04, 
0.38), and (6.71, 0.20), (top to bottom). Solid red curves are generated from 
the hydrodynamic simulation, while the dashed blue curves are generated 
by 21cmFAST. The thin solid red curve in the upper panel corresponds to a 
fully neutral universe at 2; = 9. The linear regime geometric enhancement 
of 1.87 is demarcated by the upper horizontal dotted line. 

fore the optical depth a nd STb are increased in ^ > regions 
tearkana & Loebll2005ah . 

To further explore this effect and compare our results with 
simulations, in Fig. (6] we plot the ratio of the dimensional 21 -cm 
power spectra, 6Tb{z)^ A2i{k, z), computed including peculiar ve- 
locities to those not including peculiar velocities. The thin solid red 
curve in the upper panel corresponds to a fully neutral universe at 
z = 9. The hydrodynamic simulations go down to much smaller 
scales than plotted in Fig. [5] which are more non-linear and hence 
show a larger enhancement of power, though we confirm that most 
of this is due to the evolution in the mean signal, dTb^z^. 

This enhanced 21 -cm power from non-linear peculiar veloci- 
ties obviously merits more investigation beyond the scope of this 
paper. Therefore we defer further analysis to future work. We cau- 
tion however that it is unclear how well we can estimate this en- 
hancement, due to the misuse of the 1 / (dvr /dr/H-\-l) term in eq. 
([Q. This expression assumes that dvr/dr <^ H{z) and diverges at 
dvr/dr — —H{z). To compensate for this behavior, we impose 
a maximum value of \dvr/dr\ = 0.5H(z), and confirm that our 
results are only weakly sensitive to this choice in the ^ 0.1H{z) 
- 0.7 H{z) range. Similar misuses of the mapping from real space 
to redshift s pace have already b een noted in the context of galaxy 
surveys fsee lScoccimarrdE004l and references therein). Therefore, 
if the user is interested in more accurate predictions of the 21 -cm 
signal as observed with 21 -cm interferometers, we recommend to 
turn off the velocity gradient correction in 21cmFAST, and just do 
redshift space distortions directly from the velocity field as one of 
the many necessary transformations from a comoving simulation 
box to a simulated frequency signal (e.g. Harker et al...201Qi . Mate- 
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jek et al, in preparation). The discussion below in the remainder of 
this section should not be sensitive to the above inconsistency in 
applying redshift- space corrections. 

The remaining curves in Fig.[6]do not artificially set xhi = 1, 
but use the values of xhi from the numerical simulations. Solid 
red curves are generated from the hydrodynamic simulation, while 
the dashed blue curves are generated by 21cmFAST. We confirm 
the results of MF07: that the enhancement of power due to red- 
shift space distortions vanishes in the early stages of reionization, 
and subsequently only affects small scales. Again, this is due to the 
fact that the densest regions driving most of the enhancement are 
the first ones to be covered-up by HII regions. In general, the rela- 
tive enhancement due to peculiar velocity effects is well-predicted 
by 21cmFAST on moderate to large scales, k < 1 Mpc~^, but is 
under-predicted at small scales and in a very neutral universe. This 
is attributable in part to the differences in the ionization field. Our 
ionization field algorithm, FFRT, does not prope rly capture ioniza- 
tion fronts and small-scale HII structure (Zahn et aLlboiOh . Thus 
statistics such as this one which are sensitive to small- scales are not 
accurately reproduced. 21cmFAST's under-prediction of the power 
spectrum enhancement at small scales is likely also attributable in 
part to the fact that the Nyquist frequency corresponds to larger 
scales in the 21cmFAST boxes, since these are directly computed 
on a 256^ grid, whereas the RT simulation is smoothed down from 
a 512^ grid. This means that our shot noise on small scales in higher 
than the numerical simulation's, and so fractional enhancements in 
power should be less (see below). 

Finally, we note that the curves in the bottom two panels in 
Fig. [6] dip below unity at low k. This means that during the final 
stages of reionization, peculiar velocity effects actually decrease 
power on moderate to large scales. Although not previously noted 
in the 21 -cm literature, this again can be readily understood: since 
reionization is "inside-out" on large scales, remaining neutral re- 
gions will preferentially be underdensities in the late stages. There- 
fore as the average 5 of the remaining neutral regions becomes neg- 
ative on large scales, dvr/dr becomes preferentially positive, de- 
creasing power through the 1 / {dvr / dr / H + 1) term in eq.[T] We 
confirm that both the mean signal and the large-scale power show 
this decrease due to peculiar velocities in the advanced stages of 
reionization, and it appears in both the simulations and 21cmFAST. 

2.4 Full Post-Heating Comparison of 21-cm Emission 

We now combine the terms from eq. ([T]) to provide a full compari- 
son of 21cmFAST and numerical simulations, with the assumption 
ofTs ^ Try. For the purposes of these comparisons, the two codes 
only share the same initial density field; the evolved density, ve- 
locity and ionization fields for 21cmFAST are all generated self- 
consistently as explained above. 

In Fig.|7l we plot slices through the 5Th signal, generated from 
hydrodynamic simulation, the algorithm outlined in MFO'fl and 
21cmFAST, left to right columns. Rows correspond to (z, xhi) = 
(9.00, 0.86), (7.73, 0.65), (7.04, 0.38), and (6.71, 0.20), top to bot- 
tom. 

^ This algorithm isn't exactly the same as in MF07, since partially ionized 
cells are still allowed. However this does not affect our results, since on 
scales as small as our cell size, the ionization field produced by cosmologi- 
cal RT simulations can be treated as binary, (i.e. either fully ionized or fully 
neutral; see the Appendix of Zahn et al. 2010). The ionization field becomes 
less binary in the late stages of reionization, or when hard spectra dominate 
reionization. 



As already shown in ^2.11 the density fields are modeled 
quite accurately with perturbation theory. One can also see that 
both semi-numerical schemes reproduce the large-scale HII region 
morphology (shown in black) of the RT simulations. Differences 
emerge at moderate to small scales, with the FFRT ionization al- 
gorithm of 21cmFAST generally resulting in HII regions which are 
too connected. This difference is mostly attributable to the bub- 
ble flagging algorithm; in general the "flagging-the-entire- sphere" 
algorithm of MF07 better reproduces HII morp hological st r ucture 
than the "flagging-the-central-cell" algorithm of IZahn et al.l (l2007h 
(e.g. MF07). 

In Fig.[8j we show the PDFs of 5Th for the hydrodynamic sim- 
ulation {red solid curves), 21cmFAST {blue dashed curves), and 
MF07 {magenta dotted curves). Panels correspond to {z^xni) = 
(9.00, 0.86), (7.73, 0.65), (7.04, 0.38), and (6.71, 0.20), top to bot- 
tom. The left panel was generated using the unfiltered 6Tb field 
with cell length Ax = 143/256 Mpc (effectively R - 0.35 Mpc), 
while the right panel was generated from the 6Tb field, filtered on 
i^fiiter = 5 Mpc scales. 

From the left panel, we see that we under-predict the number 
of "almost" fully-ionized cells, 6Tb < 10 mK. This can be traced 
to our algorithm for determining the partial ionized fraction in the 
remaining neutral cells. Our algorithm assumes that cells are par- 
tially ionized by sub-grid sources chewing away at their host cell's 
HI. Instead, partially ionized cells on these small scales generally 
correspond to un resolved ionizatio n fronts from non-local sources 
(see Appendix in lZahn et^ This discrepancy decreases as 

the cell size increases, since then the fraction of cells which are 
ionized by sources internal to the cell increases, and the assump- 
tion implicit in our FFRT algorithm becomes increasingly accurate. 
Aside from this, the distributions in the left panel agree very well. 
This should not be surprising, since for comparison the ionization 
efficiency of the semi-numerical schemes was chosen so that the 
mean ionized fraction at these epochs agrees with the numerical 
simulation (i.e. the spikes at 6Tb = mK matchj^- The remainder 
of the signal at ^Tb > 10 mK merely reflects the density distribu- 
tion of the neutral cells, and we have already demonstrated in ^2.11 
that our density fields match the hydrodynamic simulation quite 
welO 

The right panel of Fig. [8] shows the 6Tb distributions, 
smoothed on i^fiiter = 5 Mpc scales. As evidenced by the smaller 
relative spike at ^Tb = mK, the ionization fields on these scales 
are not as binary (i.e. either fully ionized or fully neutral) as those 
in the left panel. Thus the PDFs encode more information on the 
ionization algorithms. The top panel at xhi = 0.86, shows that 
the predicted distributions of 6Tb agree well around the mean sig- 
nal, but the semi-numerical schemes diverge from the RT in the 
wings. As mentioned before, the "flagging-the-entire- sphere" ion- 
ization algorithm from MF07 results in less connected HII regions, 
and so there are more isolated 5.0 Mpc neutral patches. This re- 
sults in an increased number of high-^Tb regions. The converse is 
true for the FFRT ionization scheme which is the fiducial setting 
of 21cmFAST. The agreement between the schemes improves as 
reionization progresses. 

Strictly speaking, the partial ionized cells do impact this calibration, but 
since most cells are either fully ionized or fully neutral, this is a small effect. 

The 6Tb distributions include the additional check of sampling the den- 
sity field of the remaining neutral cells, instead of the entire density field 
studied in ^2.11 Therefore the close agreement of the STb ^10 PDFs 
contains an additional confirmation of the accuracy of our ionization al- 
gorithms. 
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Figure 7. 6Tjj maps. The slices are generated from the hydrodynamic simulation, DexM (MF07), and 21cmFAST, left to right columns. All slices are 
Mpc on a side and 0.56 Mpc thick, and correspond to (z, xri) = (9.00, 0.86), (7.73, 0.65), (7.04, 0.38), and (6.71, 0.20), top to bottom. 
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unfiltered, Ax = 0.56 Mpc 




Figure 8. PDFs of STiy created using eq. ([T) for the hydrodynamic simulation {red solid curves), 21cmFAST {blue dashed curves), and MF07 {magenta 
dotted curves). Panels correspond to (z, xhi) = (9.00, 0.86), (7.73, 0.65), (7.04, 0.38), and (6.71, 0.20), top to bottom. The left panel was generated using 
the unfiltered ^T^ field with cell length Ax = 143/256 Mpc (effectively ~ 0.35 Mpc), while the right panel was generated from the ^T^ field, filtered on 
^filter = 5 Mpc scales. All fields invoke the simplifying assumption of Ts ^ Tj. 
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of percent leveQ At moderate to large scales, agreement is best, 
with MF07 performing slightly better than the FFRT algorithm 
which is default in 2 1 cmFASTm On smaller- scales, MF07 predicts 
too much power, while 21cmFAST under-predicts the power. It was 
shown in Zahn et al. ( 2010) that the FFRT ionization algorithm 
used in 21cmFAST over-predicts the correlation of the ionization 
and density fields on small scales, due to the fact that it operates 
directly on the evolved density field. This strong cross-correlation 
results in an under-prediction of 21 -cm power on these scales. The 
converse is true of the MF07 scheme, which although using discrete 
source halos, paints entire filtered regions as ionized, thus under- 
predicting the cross-correlation of the ionization and density fields. 
The optimal configuration for accurately estimating the 21 -cm sig- 
nal sem i-numerically is the FFRT-S scheme discussed in Zahn et alj 
boioh . set as default in the publicly- available DexIvF^. 

Most importantly, the model uncertainties of the semi- 
numerical schemes are smaller than the evolution due to reioniza- 
tion over a range Axhi 0.2. Naively therefore, one can predict 
that the semi-numerical schemes are accurate enough to estimate 
xhi from the power spectra to di < 0.1, or even better if the be- 
havior of the models are understood. However, there are many as- 
trophysical uncertainties associated with prescriptions for sources 
and sinks of ionizing photons during the epoch of reionization, and 



Figure 9. Comparison of 21 -cm power spectra obtained from the hydrody- 
namic cosmological simulation {solid red curves), and the semi-numerical 
algorithms in MF07 {dotted magenta curves) and 21cmFAST {dashed blue 
curves). Sets of curves correspond to (z, xui) = (9.00, 0.86), (7.73, 0.65), 
(7.04, 0.38), and (6.71, 0.20), top to bottom at high k. 



In Fig. |9] we compare the power spectra of these 6Tb boxes. 
Again, the hydrodynamic simulation is shown with red solid 
curves, MF07 is shown with dotted magenta curves, and 21cm- 
FAST is shown with dashed blue curves. 

At all scales, the power spectra agree with each other at the 10s 



The semi-numerical simulations show an increase in power approach- 
ing the Nyquist frequency, which is most likely just shot noise of our 256^ 
boxes. The numerical simulations do not show the same upturn, since they 
were generated on higher resolution grids (512^ for the RT and 768^ for 
the density), and subsequently smoothed down to 256^; numerical simu- 
lations generated directly on the same scale 256^ grids show simila r shot 
noise upturns in power on these scales (see Fig. 7 in lZahn et"Sll201Ql) . 

Note that the FFRT results shown here are not precisely analogous to 
those in Zahn et al. ( 2010), since there the evolved density field was taken 
from an N-body simulation, where in 21cmFAST, we self-consistently gen- 
erate the density field according to ^2.11 

http://www.astro.princeton.edu/ mesinger/DexM.html 
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it will likely be these which regulate the achievable constraints on 
xhi- Therefore it is imperative for models to be fast and be able 
to span large regions of parameter space. A single 21cmFAST re- 
alization of the 5Th fields shown in this section (generated from 
1536^ ICs) takes ^ 30 minutes to compute on a single-processor 
computer. 



3 THE SPIN TEMPERATURE 

We now relax the requirement in 321of ^ T^, and derive the 
full 21 -cm brightness temperature offset from eq. ([T]), including the 
spin temperature field. As mentioned previously, models predict 
that the heating e poch concluded well before the bulk of reioniza- 
tion, at z > 10 dFurlanettdl2006l:[c'hen & Miralda-Escudell2008l : 



ISantos et al.l l2008: 'Bac k et al.ll200 9V However, the second gener- 
ation 21 -cm interferometers, such as SKA, might be able to peak 
into this high-redshift regime of the dark ages. Furthermore, the 
astrophysical quantities at high- 2; are uncertain, and we do not re- 
ally know how robust is the assumption of Ts ^ even during 
the early stages of reionization. Therefore, for many applications, 
especially parameter studies, it is important to compute the spin 
temperature field. Unfortunately, there is currently no numerical 
simulation which includes the computationally expensive radiative 
transfer of both X-rays and Lya photons from atomically or molec- 
ularly cooled sources required to compute Ts numerically (though 
see the recent work of B ack et al. 2010, who perform RT simula- 
tions on a small subset of sources, with M > lO^^Mo). Therefore 
we cannot directly compare our the spin temperature fields to nu- 
merical simulations. 

Our derivations i n this secti o n are similar to other semi- 
analytic models (Furlanettd l2006l : IPritchard & Furlanettd l2007l : 



Santos et al. 



Santos et al. 



However, unlike ISantos et al J (120081) and 
we do not explicitly resolve the halo field as 
an intermediary step. Instead we operate directly on the evolved 
density fields, using excursion set formalism to estimate the mean 
number of sources inside spherical shells corresponding to some 
higher redshift. As discussed above, bypassing the halo field allows 
the code to be faster, with modest memory requirements. Below we 
go through our formalism in detail. 

^The spin temperature can be written as (e.g. iFurlanetto et al.l 

l2006h : 



Ty ~h XolTcX XqTj^ 
1 ~h X(y ~\~ Xc 



(5) 



where Tk is the kinetic temperature of the gas, and is the color 
temperature, which is closely coupled to the kinetic gas tempera- 
ture, Ta ~ Tk (lFieldlll959l) . There are two coupling coefficients 
in the above equation. The collisional coupling coefficient can be 
written as: 



Field dWouthuvsenll 1952l : lFieldl 19581 : WF) coupling coefficient can 
be written as: 



1.7 X 10^^(1 + z)"^5'c.Jc. 



(7) 



where Sa is a correction factor of order unity involving detailed 
atomic physics, and is the Lyman a background flux in units 
of pcn i ~^ s~ ^ Hz~^ sr~^. We compute T^ and according to 
lHiratal ( l200d) . 

According to the above equations, there are two main fields 
governing the spin temperature: (1) the kinetic temperature of the 
gas, Tk(x, z), and (2) the Lya background, Jc^(x, z). We address 
these in ^3. H and ^3.21 respectively. 



3.1 The Kinetic Temperature 

3.1.1 Evolution Equations 

To calculate the kinetic temperature, one must keep track of the in- 
homogeneous heating history of the gas. We begin by writing down 
the evolution equation for Tk (x, z) and the local ionized fraction in 
the "neutral" (i.e. outside of the ionized regions discussed in ^ I2.2D 
IGM, Xe(x,z): 



dXeip^, z') 



d^ 



dt 
dz 



;[A. 



aACXeUhfYl] 



(8) 



dt 



Sksi'^ + Xe) dz 



2Tk dm Tk dxe 
Sub dz' 1 + Xe dz' 



(9) 



where Uh = rib,o(l + + ^ni(x, z^)] is the total (H + 

He) baryonic number density at (x, z^), ep(x, z^) is the heating 
rate per baryoi^^ for process p in erg s~^, Aion is the ioniza- 
tion rate per baryon, ax is the case-A recombination coefficient, 
C = (n^) I (n)^ is the clumping factor on the scale of the simula- 
tion cell, /cs is the Boltzmann constant, /h is the hydrogen number 
fractioiH and we distinguish between the output redshift of in- 
terest, z, and some arbitrary redshift higher redshift. We also 
make the accurate assumption that single ionized helium and hy- 
drogen are ionized to the same degree, Xe(x, z), inside the mainly 
neutral IGM (e.g. Wyithe & Loeb 2003). 

In order to speed-up our calculation, we extrapolate the 
cell's density to higher redshifts assuming linear evolution from 
z (the desired output redshift at which we compute the non- 
linear density field with perturbation theory): ^ni(x, 2;^) ^ 
^ni(x, z)T)(z^)/T)(z), where T)(z) is the linear growth factor. 
In principle, we should follow the non-linear redshift evolution 
of each cell's density, ^ni(x, z^). However, linearly extrapolating 
backwards from z is a good approximation, considering that the 



0.0628 K 

AioT^ 



(6) 

where Aio = 2.85 x 10 s ^ is the spontaneous 
emission coefficient, nni, rte, and rip are the number den- 
sity of neutral hydrogen, free electrons, and protons respec- 
tively, 3id_j^^ (Tk)^ k \%[Ty^. and /^?^n(T K ) are_ taken 
^1 



from IZy gelmanI teOOSi) . IFurlanetto & Furlanettd ([20o7), and 



IFurlanetto & Furlanett respectively. The Wouthuysen 



Note that our nota t ion is different than that in iFurlanettol ( |2006|) and 
IPritchard & Furlanettol ( |2007|) , who present heating and ionization rates per 
proper volume. 

Equation ([8) ignores Helium recombinations, which is a good approx- 
imation given that most He recombining photons will cause ionizations of 
HI or Hel. 

For clarity of presentation, we will only explicitly show dependent vari- 
ables of functions on the left hand side of equations. Where is is obvious, 
we also do not explicitly show dependences on (x, z'). 
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majority of cells sized for cosmological simulations are in the lin- 
ear or quasi-linear regime at very high redshifts where heating is 
important (see, e.g., Fig. [J]) . Additionally, the evolution of struc- 
ture, and thus also of heating, is extremely rapid during this epoch, 
and is dominated by redshifts not much higher than z. Therefore, 
we can rewrite eq. ^ and eq. ^ as: 



dxe{'x.,z) dt 



dz 
dt 
dz 



dz' 



Aio 



yaACxlfHnb,o{l + zf[l + ^ni(x, 



D{z) J 



(10) 



c?Tk(x, z') 
dz' 



dt 



3/cb(1 + Xe) dz 



2Tk 



2Tk 



dD{z')/dz' 



dXe 



l + z' 3 D{z)/5r,i{yi,z)^D[z') l + Xedz'' 

(11) 

The first term in eq. (HTI) corresponds to the heat input, which for 
our purposes is dominated by the heating processes discussed be- 
low. The second term corresponds to the Hubble expansion, the 
third corresponds to adiabatic heating and cooling from structure 
formation, and the fourth corresponds to the change in the total 
number of gas particles due to ionizations. 



3.1.2 Heating and Ionization Rates 

At very high redshifts, Compton scattering between the CMB pho- 
tons and the residual free electrons sets Tk = T^. After de- 
coupling, the gas temperature evolves through adiabatic cooling, 
Ty.{z) = 2.73 X (1 + Zdec)[(l + z')/{l + Zdec)]^ whcrc the 
decoupling re dshift is given by Zdec ~ 137(nb/^V0-022)°-^ - 1 
(|Peeblesl[l993) . However, as this is o nly approximate, w e shall use 
the actual Compton heating rate (e.g. lSeager et al.l200Qh in the sum 
in eq. ([TT1) : 



3A:s(l + Xe 



1 + /He + Xe SnieC 



(T^-Tk), (12) 



where fue is the helium number fraction, is the energy den- 
sity of the CMB, and ax is the Thomson cross-section. The initial 
conditions for Xe and Tk in 21cmFAST can either be provided by 
the user, or taken from the publicly available code RECFASTrn 
(ISeageretal.lll999h . 

We now proceed to outline our procedure for estimating X- 
ray heating, which is the dominant heating process in this epoclF^. 
We compute the heating rate per particle by summing contributions 
from sources located in concentric spherical shells around (x, z'). 
First, we take the standard ansatz in assuming that sources emit 
photons with a rate proportional to the growth of the total mass frac- 
tion inside DM halos, /coii- Note that while this may not be strictly 



^ ^ http ://w w w. astro . ubc . ca/people/scott/recf ast . html 

There are two other processes that may be important. The first is shock 
heating in the IGM: with such a small temperature, a population of weak 
shocks could substantially change the thermal history. Th is appears to be 
the case in at le ast one simulatio n iGnedin & Shave j[2004[ ). although other 
simu lations ( Kuhlen et al. 2006) and analytic models ( Furlanetto & Loeb 
12004) predict much smaller effects. The second is any exotic process, such 
as dark matter decay or annihilation, that produces X-rays or hot electrons 
(■Furlanetto et al.,2006) . 



true, we are averaging over large volumes and many sources, and so 
it is probably a decent assumption in practice. Thus the total X-ray 
emission rate (in s~^) per redshift interval from luminous sources 
located between z" and z" + dz" (where z" > z') can be written 
as: 



dN^ 
dz" 



Cx/*n,peHt,0(l + ^nl 



(13) 



where the efficiency, Cx, is the number of photons per solar mass 
in stars and the remaining terms on the RHS correspond to the total 
star formation rate inside the spherical shell demarcated by z" and 
z" + dz" . Specifically, /* is the fraction of baryons converted to 
stars, and dV{z") is the comoving volume element at z" (i.e. vol- 
ume of the shell). The collapse d fraction is computed accordin g to 
a hybrid prescription, similar to lBarkana & Loebl (|2004 [2008b : 



fco\\{y^,z" ,R" 



/sT 

/ps,n: 



erfc 



y^2[Srnir.-S^"] 



(14) 

where R" is the comoving, null-geodesic distance between z' and 
z" , Smin and are the mass variances on the scales corre- 
sponding to the smallest mass sources and R", respectively, S^i is 
the evolved densit^ smoothed on scale R", which we again lin- 
early extrapolate from z: (x, z'') = 6^1' {x,z)D{z")/D{z), 
Sc is the critical linear density corresponding to virialization, 
and fsTjz" , Smin) is the mean ST collapsed fraction (with the 
I Jenkins et al.l 1200 1 normalization) while fps,n\{z" , Smin, R^^) is 
the mean PS collapse fraction, averaged over the evolved den- 
sity field, 6^1 . Therefore the normalization factor, /sT//ps,ni, 
assures that the mean collapse fraction matches th e ST collapse 
fraction (in agreement with N-body simulations ; e.g. I Jenkins et al.l 
I2OOII : iMcOuinn et alJl2007l : iTrac &Cenll2007h . while the fluctu- 
ations are sourced by the conditional PS model applied on the 
evolved density fielco- There is an implicit assumption in eq. iT3[ 
that dD{z")/dz" <^ dfcow/dz", which is quite accurate for all 
pertinent epochs, as the large-scale density evolves much slower 
than the exponential growth of the collapsed fraction. 

We assume that the X-ray luminosity of sources can be char- 
acterized with a power law of the form, Le oc (z//zyo)~°^, with 
being the lowest X-ray frequency escaping into the IGM. We can 
then write the arrival rate [i.e. number of photons s~^ Hz~^ seen 
at (x, z^)] of X-ray photons with frequency ly, from sources within 
z" and z" + dz" as: 



Although the co llapse criter ion for conditional PS was derived using the 
linear density field, Zahn et al. ( 2010) find that using the evolved density in 
eq. ( fm when computing ionization fields yields a better match to radia- 
tive transfer simulations. However, this distinction is less important in the 
context of computing X-ray flux fields than ionizing flux fields, since the 
density is more linear on the most pertinent scales R" , due to the larger 
mean free path of X-ray photons. 

Also note that eq. ( flU impHcitly assumes that the sources are evenly 
distributed within R" , which is the same assumption inherent in the ioniza- 
tion algorithm. Although the rate of change of the collapse fraction in each 
shell can be computed from PS, we hesitate to apply this prescription to 
X-rays, since the ionization algorithm which uses eq. ([14) has already been 
rigorously tested and found to agree well with RT simulations of reioniza- 
tion cZahn et al...20ia) . 
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c(x, z/, z , z") _ dNx 



dz" 



\^z" 
1 + z' 



(15) 

where the last term accounts for IGM attenuation. For computation 
efficiency, we compute the IGM X-ray mean free path through the 
meai^lGM: 



v,z ,z ) = j 



dz 



cdt 
dz 



xm{z)n{z)a{z\ z>) 



(16) 



where the photo-ionization cross-section is weighted over 

SpecieS,a(z',Z>) = /H(l-Xe)crH+/He(l-Xe)aHeI + /He^eCrHeII 

and is evaluated atz> = ^{1 + z) / {1 + z').ln practice, the contri- 
bution of Hell can be ignored. We also remind the reader that xhi 
is the volume filling factor of neutral regions during reionization. 

Finally, one obtains the X-ray heating rate per baryon, ex from 
eq. ([TTJ by integrating over frequency and z"\ 



ex(x, 2;^) ^Cy^acvQ V*^bPcrit,o(l + z'Y 

diy I — j Xl^^^ ~ El^)fhea,tfiXiai 



I 

«/ z' 



dz [l + Z ) ^ (1 + ^nl 



(18) 



Now the integrand in both integrals only depends on a single vari- 
able, and the entire frequency integral can be treated as a function 

ofz'\ 

Analogously, we can also express the ionization rate per par- 
ticle in eq. JTot as: 



Aion ( X, Z 



roo 

= / dv^^fiXiCFiFi 

Max[zyo,z^T = l] i 



i: 



dz' 



, d(j):^/dz" 



(19) 



ex(x,z^)= / diy^^{hiy - El^)fheatfiXi 



1: 



dz 



,, d (j)y^/dz" 



(17) 



where Vp is the proper, null-geodesic separation of z' and z" , and 
the frequency integral includes a sum over species, i — HI, Hel, or 
Hell, in which /^ is the number fraction, x^ is the cell's species' 
ionization fraction [which for HI and Hel is (1 — Xe), and for Hell 
is Xe], (Ji the ionization cross-section, and Ef^ is the ionization 
threshold energy of species i. The factor /heat [hf — Ef^^ Xe (x, z^)] 
is defined as the fraction of the electron energy, hh>—Ef^ , depos ited 
as heat. We use the new results of IFurlanetto & Stoevej (l2009h , to 
compute /heat, as well as /ion and /ly^ below. These fractions 
take into account the cell's local ionization state, Xe.1 as opposed to 
the global IGM value, Xe, used to calculate the optical depth in eq. 

Unfortunately, the double integral eq. ( [TtI ) is slow to evaluate 
due to the attenuation term, which depends on both redshift and fre- 
quency (eq.[T6]). To speed-up computation, we make the additional 
approximation that all photons with optical depth, rx < 1 are ab- 
sorbed, while no photons with optical depth rx > 1 are absorbed. 
Such a step-function attenuation has been shown to yield fairly ac- 
curate ionizing pho ton flux probability distributions (see Fig. 2 in 
iMesinger & Furlan etto 2009, although comparisons were limited 
to a single set of parameters). This approximation allows us to sep- 
arate the frequency and redshift integrals in eq. ([17]), removing the 
exponential attenuation term from d(j)x/dz" and setting the lower 
bound of the frequency integral to either or the frequency corre- 



sponding to an optical depth of unity, z^T=i(xe 



), whichever 



is larger. Expanding and grouping the terms, we obtain: 



\hv-E' 



/ion, HI /ion, Hel /ion, Hell i 



Fth 
^HI 



Fth 
^Hel 



^Hell 



where /ion,j [hv — Ef"^ Xe (x, z^) , j] is now the fraction of the elec- 
tron's energy going into secondary ionizations of species j, with 
the unity term inside the sum accounting for the primary ionization 
of species i. 

3.2 The Lya background 

The Lyman a background has two main contributors: X-ray exci- 
tation of HI, Ja,x; and direct stellar emission of photons between 
Lya and the Lyman limit, Ja,*. The former can easily be related 
to the X-ray heating rate, assuming that the X-ray energy injec- 
tion rate is balanced by photo ns redshifting out of Lya resonance 
(jPritchard & Furlanettoll2007h: 



Jc.,x(x, z) 



CUb 



I 



47tH{z 



Max[zyQ ,zyT- 



\l^c, Jz' 

duJ2(h^ - Ef')^-^Ux,a, , (20) 



d(l)x/dz" 



where fhya[hiy — ^*^,Xe(x, z)] is the fraction of the electron's 
energy going into Lya photons. 

Because of the high resonant optical depth of neutral hydro- 
gen, photons redshifting into any Lyman-n resonance at (x, z) 
will be absorbed in the IGM. They then quickly and locally cas- 
cade with a fraction /recycle (n) passing through Lya and inducing 
strong coupling (Hirata 2006; Pritchard & Furlanetto 2006). There- 
fore, the direct stellar emission component of the Lya background 
(in pcm~^ s~^ Hz~^ sr~^) can be estim ated with a sum ov er the 
Lyman resonance backgrounds (e.g. Barkana & Loebll2005bh : 



Note that our framework for computing the spin temperature starts to 
break down in the advanced stages of reionization. When HIT regions be- 
come sizable, there will be large sightline-to-sightline fluctuations in the 
X-ray optical depth, which can only be taken into account with approxi- 
mate radiative transfer. In our fiducial model (see below), heating has al- 
ready saturated at xm > 0.99, so this is not a concern. However, some ex- 
treme models might be able to push the heating regime well into the bulk of 
reionization; we caution the user against over-interpreting the results from 
21cmFAST in that regime. 



Ja,*(x, z)= ^ Jc(n,x, 2;) 

n — 2 

(21) 

, where the emissivity per unit redshift (no. of photons s~^ Hz~^) 
is calculated analogously to the X-ray luminosity above: 
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dt 



(22) 



Here ^(zy) is the number of photons produced per Hz per stellar 
baryon, and is evaluated at the emitted (rest frame) frequency: 



(23) 



The upper limit of the redshift integral in eq. (1211 corresponds to 
the redshift of the next Lyman resonance: 



1 + ^max(n) = (1 + z) 



(n + 1) 



(24) 



Following iBarkana & Loebl (l2005bh . we truncate the sum at 
rimax = 23, and use their Population II and Population III spectral 
models for For computational efficiency, one can rearrange 
the terms in eq. (1211 , placing the sum over Lyman transitions inside 
the redshift integral. Substituting in eq. (|22] ) and simplifying, we 
obtain: 



Ja,* (x, z) — 



47r 



poo 



-i , /\3/-, , -rR'\df coll 
1+Z) (1 + ^nl 



dz' 



/recycle (n)s(zy^) , (25) 



where the contribution from the sum over the Lyman transitions is 
a function of z\ and is zero at z' > Zmax(n = 2). 

The total Lyman a background is then just the sum of the 
above components: 



Ja,tot(x, Z) = Jc.,x(x, Z) + Ja,*(x, Z) 



(26) 



In our fiducial model, we do not explicitly take into account other 
soft-UV sources of hya such as quasars, assuming that these 
are sub-dominant to the stellar emission. However, our frame- 
work makes it simple to add additional source terms to the inte- 
grand of eq. (I2II , if the user wishes to explore such scenarios (e.g. 
l\bionteri & Gnedin.2009.) . 

3.3 Results: complete STb evolution 

All of the results in this section are from an L = 1 Gpc simula- 
tion, whose ICs are sampled on a 1800^ grid, with the final low- 
resolution boxes being 300^ (3.33 Mpc cells). Our fiducial model 
below assumes/* = 0.1, = 10^^ Mq^ (^ 1 X-ray photon per 
stellar baryonO huo = 200 eV, a = 1.5, Tvir,min = 10^ Kfor all 
sources (X-ray, Lyman a and ionizing), C = 2, i?max = 30 Mpc, 

Si 



31.d I and the stellar emissivity, s, of Pop II stars from 



iBarkana & Loebl (l2005b ) normalized to 4400 ionizing photons per 
stellar baryon. The free parameters pertaining t o the spin temper - 
ature evolution were chosen to match those in Furlanettd (I2OO6I) 
and IPritchard & Furlanettd (l2007h . to facilitate comparison. It is 
trivial to customize the code to add, for example, redshift or halo 



This number was chosen to match the total X-ray lum inosity per unit 
star formation rate at low redshifts (see iFurlanettol I2OO6I and references 
therein for details). 

This emissivity was chosen so that the midpoint of reionization is 2; ~ 
10 and the end is 2; ~ 7. 
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Figure 10. Evolution of the mean temperatures from 21cmFAST in our 
fiducial model. Solid, dashed and dotted curves show Ts , Tk and , re- 
spectively. 

mass dependences to these free parameters. The impressive length 
of the above list of uncertain astrophysical parameters (which itself 
is only a simplified description of the involved processes) serves 
well to underscore the need for a fast, portable code, capable of 
quickly scrolling through parameter space. 

We also note that the Ts calculations outlined in ^ are the 
slowest part of the 21cmFAST code (as they involve tracking evo- 
lution down to the desired redshift), and therefore should only be 
used in the regime where they are important (z > 17 in our fiducial 
model). For example, generating a 6Tb box, assuming Ts ^ T^, 
on a 300^ grid takes only a few minutes on single processor (de- 
pending on the choice of higher resolution for sampling the ICs). 
However, including the spin temperature field takes an additional 
day of computing time. Nevertheless, once the spin temperature 
evolution is computed for a given realization at all of the inter- 
mediate outputs 2Li z' > z can be used to compute 6Tb at those 
redshifts at no additional computation cost. 

Before showing detailed results, it would be useful 
to summarize the vario u s evo lutionary stages (c.f. §3.1 
in IPritchard & Furlanettd l2007h . The reader is encour- 
aged to reference the evolution of the mean temperatures 
shown in Fig. [TO] and/or view the full movie available at 
http://www.astro.princeton .edu/^mesinger/Movies/delT.mov| 
while reading below. 

(i) Collisional coupling; Tk = Ts < T^: At high redshifts, 
the IGM is dense, so the spin temperature is collisionally coupled 
to the gas kinetic temperature. The gas temperature is originally 
coupled to the CMB, but after decoupling cools adiabatically as oc 
(1 + z) faster than the CMB. The 21 -cm brightness temperature 
offset from the CMB in this regime starts at zero, when all three 
temperatures are equal, and then becomes increasingly negative as 
Ts and Tk diverge more and more from Tj. The fluctuations in 6Tb 
are driven by the density field, as collisional coupling is efficient 
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everywhere. In our fiducial model, this epoch corresponds to 100 < 
z. 

(ii) Collisional decoupling; Tr < Ts < T^: The IGM be- 
comes less dense as the Universe expands. The spin temperature 
starts to decouple from the kinetic temperature, and begins to ap- 
proach the CMB temperature again, thus 5Th starts rising towards 
zero. Decoupling from Tk occurs as a function of the local gas den- 
sity, with underdense regions decoupling first. The power spectrum 
initially steepens, as small-scale density fluctuations drive the ad- 
ditional fluctuations of the collisional coupling coefficient. As the 
spin temperature in even the overdense regions finally decouples 
from the kinetic temperature, the power spectrum flattens again, 
and the mean signal drops as Ts 0. In our fiducial model, this 
epoch corresponds to 35 < z < 100. 

(iii) Collisional decoupling WF coupling transition; Tk < 
Ts ^ Tj-. As the spin temperature throughout the IGM decouples 
from the kinetic temperature, the mean signal is faint and might 
disappear, if the first sources wait long enough to ignite. In our 
fiducial model, this transition regime doesn't really exist. In fact 
our first sources turn on before the spin temperature fully decouples 
from the kinetic temperature. 

(iv) WF coupling; Tk < Ts < T^: The first astrophysical 
sources turn on, and begin coupling the spin temperature of the 
nearby IGM to the kinetic temperature through the WF effect (Lya 
coupling). As the requirements for Lya coupling are more mod- 
est than those to heat the gas through X-ray heating, the kinetic 
temperature keeps decreasing in this epoch. The mean brightness 
temperature offset from the CMB starts becoming more negativf^ 
again and can even reach values of ^Tb < — 100 mK. In our fiducial 
model, this epoch corresponds to 25 <z < 35. 

(v) WF coupling X-ray heating transition; Tk ~ Ts < 
T^ : Lya coupling begins to saturate as most of the IGM has a spin 
temperature which is strongly coupled to the kinetic temperature. 
The mean spin temperature reaches a minimum value, and then be- 
gins increasing. A few underdense voids are left only weakly cou- 
pled as X-rays from the first sources begin heating the surrounding 
gas in earnest, raising its kinetic temperature. The 21 -cm power 
spectrum steepens dramatically as small-scale overdensities now 
host hot gas, while on large scales the gas is uniformly cold as Lya 
coupling saturates. As inhomogeneous X-ray heating continues, the 
large-scale power comes back up. In our fiducial model, this tran- 
sition occurs around z ^ 25. 

(vi) X-ray heating; Tk = Ts < T^ : X-rays start permeating 
the IGM. The fluctuations in 6Tb are now at their maximum, as 
regions close to X-ray sources are heated above the CMB temper- 
ature, STb > 0, while regions far away from sources are still very 
cold, STb < 0. A "shoulder" in the p ower spectrum, simila r to that 
seen in the epoch of reionization fe.g. lMcOuinn et al moves 
from small scales to large scales. X-rays eventually heat the entire 
IGM, and 21 -cm can only be seen in emission. The power spec- 
trum falls as this process nears completion. In our fiducial mode, 
this epoch corresponds to 18 <z <25. 

(vii) X-ray heating reionization transition; Tk = Ts > 
T^: X-rays have heated all of the IGM to temperatures above the 
CMB. The 21 -cm signal becomes insensitive to the spin temper- 
ature. Emission in 21 -cm is now at its strongest before reioniza- 
tion begins in earnest. The 21 -cm power spectrum is driven by the 



^ Note that we discuss global trends here. Locally around each Lyo; 
source, there are partially ionized regions hosting hotter gas (e.g. Cen 2006.,) . 
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Figure 11. Top panel: Evolution of the mean temperatures from 21cmFAST 
in our fiducial, = 10^^ Mq"*^ model (blue curves), and a model with a 
hundred times weaker X-ray heating, (^x = 10^^ Mq"*^ (red curves). Solid, 
dashed and dotted curves show Ts , Tk and T^ , respectively. Bottom panel: 
The corresponding evolution in Xe and xhi- 

fluctuations in the density field. In our fiducial model, this epoch 
corresponds tol6<z<18. 

(viii) Reionization: Ionizing photons from early generations of 
sources begin permeating the Universe, wiping-out the 21 -cm sig- 
nal inside ionized regions. The power spectrum initially drops on 
large scales at xhi > 0.9 as the first regions to be ionized are the 
small-scale overdensities (McOuinn et al. 2007). The mean signal 
decreases as HII regions grow, and the power spectrum is governed 
by HII morphology. This epoch can have other interesting features 
depending on the detailed evolution of the sources and sinks of 
ionizing photons, as well as feedback processes, but as the focus of 
this section is the pre-reionization regime, we shall be brief in this 
point. In our fiducial model, this epoch corresponds to 7 < z < 16. 

These milestones are fairly general, and should appear in 
most regions of astrophysical parameter space. However, the de- 
tails of the signal, as well as the precise timing and duration of 
the above epochs depends sensitively on uncertain astrophysical 
parameters. For example, note that the above epochs in our fidu- 
ci al model are shift ed t o higher redshifts than t he analogous ones 
in lFurlanett o* ('2006^ and'Pritch ard & Furlanettol (12007). This is be- 
cause the source abundances in those works were computed with 
PS, which under-predicts the abundances of Tvir > 10^ K halos 
by oy er an order of magnitude at high redshifts (e.g. iTrac & CenI 
l2007h . A similar effect is seen when com pared to the recent nu- 
merical simulations of iBaek et al. I tOld) , who were only able to 
resolve halos with mass > 10^° M©, which is 2-3 orders of mag- 
nitude away from the atomic cooling threshold at these redshifts. 
With such rare sources, they heat the gas in their \00h~^ Mpc 
boxes to above the CMB much later, at z < 9. 

Thorough investigation of the available parameter space is be- 
yond the scope of this work. However, just to briefly show an al- 
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ternate evolution, in Fig. [TT] we include a model where the X-ray 
efficiency is two orders of magnitude lower than in our fiducial 
model. As one would expect, the Lya pumping epoch is unaffected. 
However, X-ray heating is delayed by ^ 7. In such an extreme 
model, the 21 -cm signal would be seen in strong absorption against 
the CMB for a long time, and the X-ray heating epoch would over- 
lap with the early stages of the reionization epoch. 

Finally, in Fig. [TJl we show slices through our fiducial STh 
box {left), and the corresponding 3D power spectra {right). The 
slices were chosen to highlight various epochs in cosmic 21 -cm 
signal discussed above: the onset of Lya pumping, the onset of 
X-ray heating, the completion of X-ray heating, and the mid-point 
of reionization are shown from top to bottom. We encourage the 
interested reader to see more evolutionary stages through the movie 
at 'http://www.astro.princeton.edu/~mesinger/Movies/delT.mov 
When normalized to the same epoch, our power- spectrum 
evolution agrees fair l y wel l with the analytical model of 
IPritchard & Furlanettol ( 2007b. as well as its application to a 
numerical ( Santos et al. l2008h and a semi-numerical (MF07) 
simulation (Santos et al. 2009; c.f. see their Fig. 11). 

One interesting feature worth mentioning is that our reioniza- 
tion power spectra show a drop in power on large scales, which per- 
sists throughout reionization. In the advanced stages of reionization 
(^Hi < 0.9), the strongest imprint on the 21-cm power spectrum is 
from HII morphology, with a "shoulder" feature quickly propagat- 
ing from small to large sc ales, and flattening the power spectrum 
(e.g. iMcOuinn et al ] |2007l : see Fig. |9]). However there should still 
be a drop in large-scale power beyond either the HII bubble scale, 
or the photon mean free path in the ionized IGlvQ whichever is 
smallest. This is an interesting feature from which one can de- 
duce the ionizing photon mean free path in the late stages of reion - 
ization dMesinger & Diikstrall2008l : iFurlanetto & Mesingejl2009h . 
And since it occurs at A; < 0.1 Mpc~^ (for our fiducial choice of 
parameters), it is beyond the dynamic range of present-day numer- 
ical simulations hoping to resolve atomically-cooled source halos. 



4 CONCLUSIONS 

We introduce a powerful new semi-numeric modeling tool, 21cm- 
FAST, designed to efficiently simulate the cosmological 21-cm sig- 
nal. Our approach uses perturbation theory, excursion set formal- 
ism, and analytic prescriptions to generate evolved density, ion- 
ization, peculiar velocity, and spin temperature fields, which it 
then combines to compute the 21-cm brightness temperature. This 
code is based on the semi-numerical simulation, DexM, (MF07). 
However, here we bypasses the halo finder and operate directly 
on the evolved density field, thereby increasing the speed and de- 
creasing memory requirements. In the post-heating regime, 21cm- 
FAST can generate a realization in a few minutes on a single 
processor, compared to many days on > 1000-node supercom- 
puting cluster required to generate the same resolution boxes 
using state-of-the-art numerical simulations. 21cmFAST realiza- 
tions in the pre-heating regime require ~ one day of computa- 
tion time. Conversely, RT simulations of the pre-heating regime 
which resolve most sources currently do not exist, as they are 
too computationally expensive. Our code is publicly available at 
[httpV/www.astro.princeton.edu/^mesinger/Sim.html] 



■^^ Note that RT simulations of reionization generally do not explicitly in- 
clude an effective ionizing photon mean free path from unresolved LLSs. 



We compare maps, PDFs and power spectra from 21cmFAST, 
with corresponding ones from the hydrodynamic numerical simu- 
lations of Trac et al. (2008), generated from the same initial con- 
ditions. We find good agreement with the numerical simulation 
on scales pertinent to the upcoming observations (> 1 Mpc). The 
power spectra from 21cmFAST agree with those generated from the 
numerical simulation to within 10s of percent down to the Nyquist 
frequency. 

We find evidence that non-linear peculiar velocity effects en- 
hance the 21-cm power spectrum, beyond the expected geometric, 
linear value. This enhancement quickly diminishes during the onset 
of reionization, remaining only on small scales at xhi ^0-7. In- 
terestingly, we also find that the large-scale power is decreased as a 
result of peculiar velocities in the advanced stages of reionization. 
The reason for this is due to the "inside-out" nature of reionization 
on large- scales: the remaining HI regions are preferentially under- 
dense, in which peculiar velocities decrease the 21-cm optical depth 
and brightness temperature. 

Our code can also simulate the pre-reionization regime, in- 
cluding the astrophysical processes of X-ray heating and the WF 
effect. We show results from a 1 Gpc simulation which tracks the 
cosmic 21-cm signal down from z — 250, highlighting the various 
interesting epochs. 

There are several large 21-cm interferometers scheduled to be- 
come operational soon. Interpreting their upcoming data will be 
difficult since we know very little about the astrophysical processes 
at high redshifts. Furthermore, there is an enormous range of scales 
involved, making numerical simulations too slow for efficient pa- 
rameter exploration. 21cmFAST is not. 
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